############################################################################################################################
########Inclusion, recognition, and inter-group comparisons: The effects of power-sharing institutions on grievances########
#####Individual-level analyses (main and robustness checks)#####
#####Andreas Juon, 2023 (Journal of Conflict Research)#####
############################################################################################################################

####################################################################################
#########STEP 1: LIBRARIES, DATA IMPORT, DEFINITIONS#########
####################################################################################

setwd("XXXXX")#set to directory that contains the data files

#######1.1 Libraries#######
library(lme4)
library(stargazer)
library(ggplot2)
library(grid)
library(gridExtra)
library(modeest)
library(dotwhisker)
library(boot)

#######1.2 Reading Data#######
grievances_ind <- read.csv(file="./juon_jcr_grievances_data_ind.csv")
grievances_ind$cowcode <- as.factor(grievances_ind$cowcode)
grievances_ind$gwgroupid <- as.factor(grievances_ind$gwgroupid)
grievances_ind$cowcode_year <- as.factor(grievances_ind$cowcode_year)
grievances_ind$gwgroupid_year <- as.factor(grievances_ind$gwgroupid_year)
grievances_ind$region <- as.factor(grievances_ind$region)
grievances_ind$survey <- as.factor(grievances_ind$survey)

#######1.3 Defining variable lists#######
main1 <- c("state_control","included_nc")
main2 <- c("state_control","state_control * ps_corp","state_control * ps_lib")
main3 <- c("state_control","state_control * ps_corp","state_control * ps_lib","state_control * ps_corp_all_diff_neg")
main4 <- c("state_control","state_control * ps_corp","state_control * ps_lib","state_control * tek_ps_corp_diff_neg")
main2_int <- c("state_control","state_control * ps_corp","state_control * ps_lib")
main3_int <- c("state_control","state_control * ps_corp","state_control * ps_lib","state_control * ps_corp_all_diff_neg * interest_pol_d")
main4_int <- c("state_control","state_control * ps_corp","state_control * ps_lib","state_control * tek_ps_corp_diff_neg * interest_pol_d")
main1_stocks <- c("state_control","included_stock")
main2_stocks <- c("state_control","state_control * ps_corp_stock","state_control * ps_lib_stock")
re <- c("1", "(1|cowcode)")#, "(1|gwgroupid)"
re_i <- c("(1|cowcode_year)")
re_i2 <- c("(1|gwgroupid_year)")
re_i3 <- c("(1|gwgroupid)")
controls_gc <- c("size","other","ongoing_recent_contestation","cpi","politya","lgdppc","fractionalization1","region")
controls_gc_noreg <- c("size","other","ongoing_recent_contestation","cpi","politya","lgdppc","fractionalization1")
controls_gc_min <- c("size","geo_conc","ongoing_contestation","recent_contestation","other")#mminority
controls_gc2 <- c("autonomy","low_ratio","high_ratio","cleavage_mean")
controls_i <- c("age","gender","educ_high","interest_pol_d")
ac1 <- c("autonomy","low_ratio","high_ratio","cleavage_mean")
ac2 <- c("country_contestation","tek_contestation","only_one_nc")
ac3 <- c("region")

#######1.4 Defining variable reporting order for stargazer#######
vars.order <- c("state_control","included_nc","ps_corp","ps_corp_all_diff_neg","tek_ps_corp_diff_neg","state_control:ps_corp","ps_corp:state_control","state_control:ps_corp_all_diff_neg","ps_corp_all_diff_neg:state_control","state_control:tek_ps_corp_diff_neg","tek_ps_corp_diff_neg:state_control","ps_lib","state_control:ps_lib","ps_lib:state_control")
vars.order_stocks <- c("state_control","included_stock","ps_corp_stock","state_control:ps_corp_stock","ps_corp_stock:state_control","ps_lib_stock","state_control:ps_lib_stock","ps_lib_stock:state_control")


####################################################################################
#########STEP 2: MAIN MODELS#########
####################################################################################

#######2.1 Main models (main article, table 1 / appendix 2, table A3)#######
#a) unsat_gov
m1_unsat_gov <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
m2_unsat_gov <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
m3_unsat_gov <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
m4_unsat_gov <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
#b) disc_grp
m1_disc_grp <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50))
m2_disc_grp <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50))
m3_disc_grp <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50))
m4_disc_grp <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50))
stargazer(m1_unsat_gov, m2_unsat_gov, m3_unsat_gov, m4_unsat_gov, m1_disc_grp, m2_disc_grp, m3_disc_grp, m4_disc_grp, type="text", style = "ajps", order=paste0("^", vars.order , "$"), omit = c("survey","region"), dep.var.labels = c("Government dissatisfaction", "Feeling discriminated"), notes = c("* p<0.1; ** p<0.05; *** p<0.01. Survey- and region-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Core","Included non-core","Corporate PSI (group)","Corporate PSI (difference, domestic)","Corporate PSI (difference, TEK)","Core x corporate PSI (group)","Core x corporate PSI (difference, domestic)","Core x corporate PSI (difference, TEK)","Liberal PSI (group)","Core x liberal PSI (group)","Relative size","non-EPR group","Recent contestation","CPI","Polity (normalized)","GDP p.c. (logged)","Ethnic fractionalization","Age","Female","High education","Political interest"))


####################################################################################
#########STEP 3: ROBUSTNESS CHECKS#########
####################################################################################

#######3.1.1 Alternative dependent variable: Central government items (appendix 4.1.1, table A7)#######
m1_unsat_cgov <- glmer(as.formula(paste("unsat_cgov", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
m2_unsat_cgov <- glmer(as.formula(paste("unsat_cgov", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
m3_unsat_cgov <- glmer(as.formula(paste("unsat_cgov", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
m4_unsat_cgov <- glmer(as.formula(paste("unsat_cgov", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
stargazer(m1_unsat_cgov, m2_unsat_cgov, m3_unsat_cgov, m4_unsat_cgov, type="text", style = "ajps", order=paste0("^", vars.order , "$"), omit = c("survey","region"), dep.var.labels = c("Executive dissatisfaction"), notes = c("* p<0.1; ** p<0.05; *** p<0.01. Survey- and region-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Core","Included non-core","Corporate PSI (group)","Corporate PSI (difference, domestic)","Corporate PSI (difference, TEK)","Core x corporate PSI (group)","Core x corporate PSI (difference, domestic)","Core x corporate PSI (difference, TEK)","Liberal PSI (group)","Core x liberal PSI (group)","Relative size","non-EPR group","Recent contestation","CPI","Polity (normalized)","GDP p.c. (logged)","Ethnic fractionalization","Age","Female","High education","Political interest"))

#######3.1.2 Sequential exclusion of each uniquely worded question item (appendix 4.1.2, figure A6)#######

###3.1.2.1 model 3###
#a) jacknife estimation
m3_unsat_gov_jacknife <- data.frame()
survey_list <- unique(subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50)$surveywave)
for (i in survey_list){
  m3_unsat_gov_temp <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50 & surveywave!=i))
  jacknife_results_temp <- data.frame(survey = i)
  coefs <- data.frame(summary(m3_unsat_gov_temp)$coefficients)
  coefs$variable <- rownames(coefs)
  coefs_ps_corp <- subset(coefs, variable == "ps_corp")
  coefs_ps_corp_all_diff_neg <- subset(coefs, variable == "ps_corp_all_diff_neg")
  jacknife_results_temp$jacknife_ps_corp <- as.numeric(round(coefs_ps_corp$Estimate, digits=10))
  jacknife_results_temp$jacknife_ps_corp_all_diff_neg <- as.numeric(round(coefs_ps_corp_all_diff_neg$Estimate, digits=10))
  m3_unsat_gov_jacknife <- rbind(m3_unsat_gov_jacknife, jacknife_results_temp)
}
#b) plot
m3_unsat_gov_jacknife_ps_corp <- ggplot(aes(x = jacknife_ps_corp), data = m3_unsat_gov_jacknife) + geom_histogram(binwidth = 0.05) + 
  ylab("Frequency") + xlab("Coefficient estimate") + ggtitle("Corporate PSI (group) [model 3]") +
  theme_bw() + theme(legend.position = "bottom", axis.text = element_text(size=7), text=element_text(family="Times New Roman"))
m3_unsat_gov_ps_corp_all_diff_neg <- ggplot(aes(x = jacknife_ps_corp_all_diff_neg), data = m3_unsat_gov_jacknife) + geom_histogram(binwidth = 0.05) + 
  ylab("Frequency") + xlab("Coefficient estimate") + ggtitle("Corporate PSI (difference, domestic) [model 3]") +
  theme_bw() + theme(legend.position = "bottom", axis.text = element_text(size=7), text=element_text(family="Times New Roman"))
m3_jacknife_plot <- grid.arrange(m3_unsat_gov_jacknife_ps_corp+ theme(plot.title = element_text(size=10)), m3_unsat_gov_ps_corp_all_diff_neg+ theme(plot.title = element_text(size=10)), nrow=1, ncol=2)

###3.1.2.2 model 4###
#a) jacknife estimation
m4_unsat_gov_jacknife <- data.frame()
for (i in survey_list){
  m4_unsat_gov_temp <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50 & surveywave!=i))
  jacknife_results_temp <- data.frame(survey = i)
  coefs <- data.frame(summary(m4_unsat_gov_temp)$coefficients)
  coefs$variable <- rownames(coefs)
  coefs_ps_corp <- subset(coefs, variable == "ps_corp")
  coefs_tek_ps_corp_diff_neg <- subset(coefs, variable == "tek_ps_corp_diff_neg")
  jacknife_results_temp$jacknife_ps_corp <- as.numeric(round(coefs_ps_corp$Estimate, digits=10))
  jacknife_results_temp$jacknife_tek_ps_corp_diff_neg <- as.numeric(round(coefs_tek_ps_corp_diff_neg$Estimate, digits=10))
  m4_unsat_gov_jacknife <- rbind(m4_unsat_gov_jacknife, jacknife_results_temp)
}
#b) plot
m4_unsat_gov_jacknife_ps_corp <- ggplot(aes(x = jacknife_ps_corp), data = m4_unsat_gov_jacknife) + geom_histogram(binwidth = 0.05) + 
  ylab("Frequency") + xlab("Coefficient estimate") + ggtitle("Corporate PSI (group) [model 4]") +
  theme_bw() + theme(legend.position = "bottom", axis.text = element_text(size=7), text=element_text(family="Times New Roman"))
m4_unsat_gov_tek_ps_corp_diff_neg <- ggplot(aes(x = jacknife_tek_ps_corp_diff_neg), data = m4_unsat_gov_jacknife) + geom_histogram(binwidth = 0.05) + 
  ylab("Frequency") + xlab("Coefficient estimate") + ggtitle("Corporate PSI (difference, TEK) [model 4]") +
  theme_bw() + theme(legend.position = "bottom", axis.text = element_text(size=7), text=element_text(family="Times New Roman"))
m4_jacknife_plot <- grid.arrange(m4_unsat_gov_jacknife_ps_corp+ theme(plot.title = element_text(size=10)), m4_unsat_gov_tek_ps_corp_diff_neg+ theme(plot.title = element_text(size=10)), nrow=1, ncol=2)

###3.1..2.3 model 7###
#a) jacknife estimation
m3_disc_grp_jacknife <- data.frame()
survey_list <- unique(subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50)$surveywave)
for (i in survey_list){
  m3_disc_grp_temp <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50 & surveywave!=i))
  jacknife_results_temp <- data.frame(survey = i)
  coefs <- data.frame(summary(m3_disc_grp_temp)$coefficients)
  coefs$variable <- rownames(coefs)
  coefs_ps_corp <- subset(coefs, variable == "ps_corp")
  coefs_ps_corp_all_diff_neg <- subset(coefs, variable == "ps_corp_all_diff_neg")
  jacknife_results_temp$jacknife_ps_corp <- as.numeric(round(coefs_ps_corp$Estimate, digits=10))
  jacknife_results_temp$jacknife_ps_corp_all_diff_neg <- as.numeric(round(coefs_ps_corp_all_diff_neg$Estimate, digits=10))
  m3_disc_grp_jacknife <- rbind(m3_disc_grp_jacknife, jacknife_results_temp)
}
#b) plot
m3_disc_grp_jacknife_ps_corp <- ggplot(aes(x = jacknife_ps_corp), data = m3_disc_grp_jacknife) + geom_histogram(binwidth = 0.05) + 
  ylab("Frequency") + xlab("Coefficient estimate") + ggtitle("Corporate PSI (group) [model 7]") +
  theme_bw() + theme(legend.position = "bottom", axis.text = element_text(size=7), text=element_text(family="Times New Roman"))
m3_disc_grp_ps_corp_all_diff_neg <- ggplot(aes(x = jacknife_ps_corp_all_diff_neg), data = m3_disc_grp_jacknife) + geom_histogram(binwidth = 0.05) + 
  ylab("Frequency") + xlab("Coefficient estimate") + ggtitle("Corporate PSI (difference, domestic) [model 7]") +
  theme_bw() + theme(legend.position = "bottom", axis.text = element_text(size=7), text=element_text(family="Times New Roman"))
m7_jacknife_plot <- grid.arrange(m3_disc_grp_jacknife_ps_corp+ theme(plot.title = element_text(size=10)), m3_disc_grp_ps_corp_all_diff_neg+ theme(plot.title = element_text(size=10)), nrow=1, ncol=2)

###3.1.2.3 model 8###
#a) jacknife estimation#
m4_disc_grp_jacknife <- data.frame()
for (i in survey_list){
  m4_disc_grp_temp <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50 & surveywave!=i))
  jacknife_results_temp <- data.frame(survey = i)
  coefs <- data.frame(summary(m4_disc_grp_temp)$coefficients)
  coefs$variable <- rownames(coefs)
  coefs_ps_corp <- subset(coefs, variable == "ps_corp")
  coefs_tek_ps_corp_diff_neg <- subset(coefs, variable == "tek_ps_corp_diff_neg")
  jacknife_results_temp$jacknife_ps_corp <- as.numeric(round(coefs_ps_corp$Estimate, digits=10))
  jacknife_results_temp$jacknife_tek_ps_corp_diff_neg <- as.numeric(round(coefs_tek_ps_corp_diff_neg$Estimate, digits=10))
  m4_disc_grp_jacknife <- rbind(m4_disc_grp_jacknife, jacknife_results_temp)
}
#b) plot
m4_disc_grp_jacknife_ps_corp <- ggplot(aes(x = jacknife_ps_corp), data = m4_disc_grp_jacknife) + geom_histogram(binwidth = 0.05) + 
  ylab("Frequency") + xlab("Coefficient estimate") + ggtitle("Corporate PSI (group) [model 8]") +
  theme_bw() + theme(legend.position = "bottom", axis.text = element_text(size=7), text=element_text(family="Times New Roman"))
m4_disc_grp_tek_ps_corp_diff_neg <- ggplot(aes(x = jacknife_tek_ps_corp_diff_neg), data = m4_disc_grp_jacknife) + geom_histogram(binwidth = 0.05) + 
  ylab("Frequency") + xlab("Coefficient estimate") + ggtitle("Corporate PSI (difference, TEK) [model 8]") +
  theme_bw() + theme(legend.position = "bottom", axis.text = element_text(size=7), text=element_text(family="Times New Roman"))
m8_jacknife_plot <- grid.arrange(m4_disc_grp_jacknife_ps_corp+ theme(plot.title = element_text(size=10)), m4_disc_grp_tek_ps_corp_diff_neg+ theme(plot.title = element_text(size=10)), nrow=1, ncol=2)


#######3.2 Alternative thresholds of correct identification (individual)#######

###3.2.1 60% (appendix 4.2.1, table A8)###
#a) unsat_gov
m1_unsat_gov_prob0.6 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat0.6 == 1 & group_s_number_unsat0.6 >= 50))
m2_unsat_gov_prob0.6 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat0.6 == 1 & group_s_number_unsat0.6 >= 50))
m3_unsat_gov_prob0.6 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat0.6 == 1 & group_s_number_unsat0.6 >= 50))
m4_unsat_gov_prob0.6 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat0.6 == 1 & group_s_number_unsat0.6 >= 50))
#b) disc_grp
m1_disc_grp_prob0.6 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc0.6 == 1 & group_s_number_disc0.6 >= 50))
m2_disc_grp_prob0.6 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc0.6 == 1 & group_s_number_disc0.6 >= 50))
m3_disc_grp_prob0.6 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc0.6 == 1 & group_s_number_disc0.6 >= 50))
m4_disc_grp_prob0.6 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc0.6 == 1 & group_s_number_disc0.6 >= 50))
stargazer(m1_unsat_gov_prob0.6, m2_unsat_gov_prob0.6, m3_unsat_gov_prob0.6, m4_unsat_gov_prob0.6, m1_disc_grp_prob0.6, m2_disc_grp_prob0.6, m3_disc_grp_prob0.6, m4_disc_grp_prob0.6, type="text", style = "ajps", order=paste0("^", vars.order , "$"), omit = c("survey","region"), dep.var.labels = c("Government dissatisfaction", "Feeling discriminated"), notes = c("* p<0.1; ** p<0.05; *** p<0.01. Survey- and region-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Core","Included non-core","Corporate PSI (group)","Corporate PSI (difference, domestic)","Corporate PSI (difference, TEK)","Core x corporate PSI (group)","Core x corporate PSI (difference, domestic)","Core x corporate PSI (difference, TEK)","Liberal PSI (group)","Core x liberal PSI (group)","Relative size","non-EPR group","Recent contestation","CPI","Polity (normalized)","GDP p.c. (logged)","Ethnic fractionalization","Age","Female","High education","Political interest"))

###3.2.2 70% (appendix 4.2.2, table A9)
#a) unsat_gov
m1_unsat_gov_prob0.7 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat0.7 == 1 & group_s_number_unsat0.7 >= 50))
m2_unsat_gov_prob0.7 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat0.7 == 1 & group_s_number_unsat0.7 >= 50))
m3_unsat_gov_prob0.7 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat0.7 == 1 & group_s_number_unsat0.7 >= 50))
m4_unsat_gov_prob0.7 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat0.7 == 1 & group_s_number_unsat0.7 >= 50))
#b) disc_grp
m1_disc_grp_prob0.7 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc0.7 == 1 & group_s_number_disc0.7 >= 50))
m2_disc_grp_prob0.7 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc0.7 == 1 & group_s_number_disc0.7 >= 50))
m3_disc_grp_prob0.7 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc0.7 == 1 & group_s_number_disc0.7 >= 50))
m4_disc_grp_prob0.7 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc0.7 == 1 & group_s_number_disc0.7 >= 50))
stargazer(m1_unsat_gov_prob0.7, m2_unsat_gov_prob0.7, m3_unsat_gov_prob0.7, m4_unsat_gov_prob0.7, m1_disc_grp_prob0.7, m2_disc_grp_prob0.7, m3_disc_grp_prob0.7, m4_disc_grp_prob0.7, type="text", style = "ajps", order=paste0("^", vars.order , "$"), omit = c("survey","region"), dep.var.labels = c("Government dissatisfaction", "Feeling discriminated"), notes = c("* p<0.1; ** p<0.05; *** p<0.01. Survey- and region-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Core","Included non-core","Corporate PSI (group)","Corporate PSI (difference, domestic)","Corporate PSI (difference, TEK)","Core x corporate PSI (group)","Core x corporate PSI (difference, domestic)","Core x corporate PSI (difference, TEK)","Liberal PSI (group)","Core x liberal PSI (group)","Relative size","non-EPR group","Recent contestation","CPI","Polity (normalized)","GDP p.c. (logged)","Ethnic fractionalization","Age","Female","High education","Political interest"))

###3.2.3 90% (appendix 4.2.3, table A10)
#a) unsat_gov
m1_unsat_gov_prob0.9 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat0.9 == 1 & group_s_number_unsat0.9 >= 50))
m2_unsat_gov_prob0.9 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat0.9 == 1 & group_s_number_unsat0.9 >= 50))
m3_unsat_gov_prob0.9 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat0.9 == 1 & group_s_number_unsat0.9 >= 50))
m4_unsat_gov_prob0.9 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat0.9 == 1 & group_s_number_unsat0.9 >= 50))
#b) disc_grp
m1_disc_grp_prob0.9 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc0.9 == 1 & group_s_number_disc0.9 >= 50))
m2_disc_grp_prob0.9 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc0.9 == 1 & group_s_number_disc0.9 >= 50))
m3_disc_grp_prob0.9 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc0.9 == 1 & group_s_number_disc0.9 >= 50))
m4_disc_grp_prob0.9 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc0.9 == 1 & group_s_number_disc0.9 >= 50))
stargazer(m1_unsat_gov_prob0.9, m2_unsat_gov_prob0.9, m3_unsat_gov_prob0.9, m4_unsat_gov_prob0.9, m1_disc_grp_prob0.9, m2_disc_grp_prob0.9, m3_disc_grp_prob0.9, m4_disc_grp_prob0.9, type="text", style = "ajps", order=paste0("^", vars.order , "$"), omit = c("survey","region"), dep.var.labels = c("Government dissatisfaction", "Feeling discriminated"), notes = c("* p<0.1; ** p<0.05; *** p<0.01. Survey- and region-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Core","Included non-core","Corporate PSI (group)","Corporate PSI (difference, domestic)","Corporate PSI (difference, TEK)","Core x corporate PSI (group)","Core x corporate PSI (difference, domestic)","Core x corporate PSI (difference, TEK)","Liberal PSI (group)","Core x liberal PSI (group)","Relative size","non-EPR group","Recent contestation","CPI","Polity (normalized)","GDP p.c. (logged)","Ethnic fractionalization","Age","Female","High education","Political interest"))

###3.2.4 all respondents, irrespective of number (appendix 4.2.4, table A11)###
#a) unsat_gov
m1_unsat_gov_allresp <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat0 == 1))
m2_unsat_gov_allresp <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat0 == 1))
m3_unsat_gov_allresp <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat0 == 1))
m4_unsat_gov_allresp <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat0 == 1))
#b) disc_grp
m1_disc_grp_allresp <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc0 == 1))
m2_disc_grp_allresp <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc0 == 1))
m3_disc_grp_allresp <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc0 == 1))
m4_disc_grp_allresp <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc0 == 1))
stargazer(m1_unsat_gov_allresp, m2_unsat_gov_allresp, m3_unsat_gov_allresp, m4_unsat_gov_allresp, m1_disc_grp_allresp, m2_disc_grp_allresp, m3_disc_grp_allresp, m4_disc_grp_allresp, type="text", style = "ajps", order=paste0("^", vars.order , "$"), omit = c("survey","region"), dep.var.labels = c("Government dissatisfaction", "Feeling discriminated"), notes = c("* p<0.1; ** p<0.05; *** p<0.01. Survey- and region-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Core","Included non-core","Corporate PSI (group)","Corporate PSI (difference, domestic)","Corporate PSI (difference, TEK)","Core x corporate PSI (group)","Core x corporate PSI (difference, domestic)","Core x corporate PSI (difference, TEK)","Liberal PSI (group)","Core x liberal PSI (group)","Relative size","non-EPR group","Recent contestation","CPI","Polity (normalized)","GDP p.c. (logged)","Ethnic fractionalization","Age","Female","High education","Political interest"))

###3.2.5 Only non-core groups (appendix 4.2.5, table A12)###
#a) unsat_gov
m1_unsat_gov_nc <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50 & state_control == 0))
m2_unsat_gov_nc <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50 & state_control == 0))
m3_unsat_gov_nc <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50 & state_control == 0))
m4_unsat_gov_nc <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50 & state_control == 0))
#b) disc_grp
m1_disc_grp_nc <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50 & state_control == 0))
m2_disc_grp_nc <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50 & state_control == 0))
m3_disc_grp_nc <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50 & state_control == 0))
m4_disc_grp_nc <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50 & state_control == 0))
stargazer(m1_unsat_gov_nc, m2_unsat_gov_nc, m3_unsat_gov_nc, m4_unsat_gov_nc, m1_disc_grp_nc, m2_disc_grp_nc, m3_disc_grp_nc, m4_disc_grp_nc, type="text", style = "ajps", order=paste0("^", vars.order , "$"), omit = c("survey","region"), dep.var.labels = c("Government dissatisfaction", "Feeling discriminated"), notes = c("* p<0.1; ** p<0.05; *** p<0.01. Survey- and region-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Included non-core","Corporate PSI (group)","Corporate PSI (difference, domestic)","Corporate PSI (difference, TEK)","Liberal PSI (group)","Relative size","non-EPR group","Recent contestation","CPI","Polity (normalized)","GDP p.c. (logged)","Ethnic fractionalization","Age","Female","High education","Political interest"))


#######3.3 Different contexts (sample subsets)#######

###3.3.1 Only democratic contexts (appendix 4.3.1, table A13)###
#a) unsat_gov
m1_unsat_gov_dem <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50 & politya >= 0.8))
m2_unsat_gov_dem <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50 & politya >= 0.8))
m3_unsat_gov_dem <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50 & politya >= 0.8))
m4_unsat_gov_dem <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50 & politya >= 0.8))
#b) disc_grp
m1_disc_grp_dem <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50 & politya >= 0.8))
m2_disc_grp_dem <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50 & politya >= 0.8))
m3_disc_grp_dem <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50 & politya >= 0.8))
m4_disc_grp_dem <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50 & politya >= 0.8))
stargazer(m1_unsat_gov_dem, m2_unsat_gov_dem, m3_unsat_gov_dem, m4_unsat_gov_dem, m1_disc_grp_dem, m2_disc_grp_dem, m3_disc_grp_dem, m4_disc_grp_dem, type="text", style = "ajps", order=paste0("^", vars.order , "$"), omit = c("survey","region"), dep.var.labels = c("Government dissatisfaction", "Feeling discriminated"), notes = c("* p<0.1; ** p<0.05; *** p<0.01. Survey- and region-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Core","Included non-core","Corporate PSI (group)","Corporate PSI (difference, domestic)","Corporate PSI (difference, TEK)","Core x corporate PSI (group)","Core x corporate PSI (difference, domestic)","Core x corporate PSI (difference, TEK)","Liberal PSI (group)","Core x liberal PSI (group)","Relative size","non-EPR group","Recent contestation","CPI","Polity (normalized)","GDP p.c. (logged)","Ethnic fractionalization","Age","Female","High education","Political interest"))

###3.3.2 Only post-conflict contexts (appendix 4.3.2, table A14)###
#a) unsat_gov
m1_unsat_gov_postc <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50 & country_contestation == 1))
m2_unsat_gov_postc <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50 & country_contestation == 1))
m3_unsat_gov_postc <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50 & country_contestation == 1))
m4_unsat_gov_postc <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50 & country_contestation == 1))
#b) disc_grp
m1_disc_grp_postc <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50 & country_contestation == 1))
m2_disc_grp_postc <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50 & country_contestation == 1))
m3_disc_grp_postc <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50 & country_contestation == 1))
m4_disc_grp_postc <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50 & country_contestation == 1))
stargazer(m1_unsat_gov_postc, m2_unsat_gov_postc, m3_unsat_gov_postc, m4_unsat_gov_postc, m1_disc_grp_postc, m2_disc_grp_postc, m3_disc_grp_postc, m4_disc_grp_postc, type="text", style = "ajps", order=paste0("^", vars.order , "$"), omit = c("survey","region"), dep.var.labels = c("Government dissatisfaction", "Feeling discriminated"), notes = c("* p<0.1; ** p<0.05; *** p<0.01. Survey- and region-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Core","Included non-core","Corporate PSI (group)","Corporate PSI (difference, domestic)","Corporate PSI (difference, TEK)","Core x corporate PSI (group)","Core x corporate PSI (difference, domestic)","Core x corporate PSI (difference, TEK)","Liberal PSI (group)","Core x liberal PSI (group)","Relative size","non-EPR group","Recent contestation","CPI","Polity (normalized)","GDP p.c. (logged)","Ethnic fractionalization","Age","Female","High education","Political interest"))


#######3.4 Additional controls#######

###3.4.1 Other horizontal inequalities (appendix 4.4.1, table A15)###
#a) unsat_gov
m1_unsat_gov_ac1 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, ac1, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
m2_unsat_gov_ac1 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, ac1, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
m3_unsat_gov_ac1 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, ac1, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
m4_unsat_gov_ac1 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, ac1, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
#b) disc_grp
m1_disc_grp_ac1 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, ac1, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50))
m2_disc_grp_ac1 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, ac1, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50))
m3_disc_grp_ac1 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, ac1, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50))
m4_disc_grp_ac1 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, ac1, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50))
stargazer(m1_unsat_gov_ac1, m2_unsat_gov_ac1, m3_unsat_gov_ac1, m4_unsat_gov_ac1, m1_disc_grp_ac1, m2_disc_grp_ac1, m3_disc_grp_ac1, m4_disc_grp_ac1, type="text", style = "ajps", order=paste0("^", vars.order , "$"), omit = c("survey","region"), dep.var.labels = c("Government dissatisfaction", "Feeling discriminated"), notes = c("* p<0.1; ** p<0.05; *** p<0.01. Survey- and region-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Core","Included non-core","Corporate PSI (group)","Corporate PSI (difference, domestic)","Corporate PSI (difference, TEK)","Core x corporate PSI (group)","Core x corporate PSI (difference, domestic)","Core x corporate PSI (difference, TEK)","Liberal PSI (group)","Core x liberal PSI (group)","Relative size","non-EPR group","Recent contestation","CPI","Polity (normalized)","GDP p.c. (logged)","Ethnic fractionalization","Autonomy","Low ratio","High ratio","Cultural cleavages","Age","Female","High education","Political interest"))

###3.4.2 Alternative causal paths (appendix 4.4.2, table A16)###
#a) unsat_gov
m1_unsat_gov_ac2 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, ac2, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
m2_unsat_gov_ac2 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, ac2, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
m3_unsat_gov_ac2 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, ac2, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
m4_unsat_gov_ac2 <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, ac2, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
#b) disc_grp
m1_disc_grp_ac2 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main1,controls_gc, ac2, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50))
m2_disc_grp_ac2 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main2,controls_gc, ac2, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50))
m3_disc_grp_ac2 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main3,controls_gc, ac2, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50))
m4_disc_grp_ac2 <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main4,controls_gc, ac2, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50))
stargazer(m1_unsat_gov_ac2, m2_unsat_gov_ac2, m3_unsat_gov_ac2, m4_unsat_gov_ac2, m1_disc_grp_ac2, m2_disc_grp_ac2, m3_disc_grp_ac2, m4_disc_grp_ac2, type="text", style = "ajps", order=paste0("^", vars.order , "$"), omit = c("survey","region"), dep.var.labels = c("Government dissatisfaction", "Feeling discriminated"), notes = c("* p<0.1; ** p<0.05; *** p<0.01. Survey- and region-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Core","Included non-core","Corporate PSI (group)","Corporate PSI (difference, domestic)","Corporate PSI (difference, TEK)","Core x corporate PSI (group)","Core x corporate PSI (difference, domestic)","Core x corporate PSI (difference, TEK)","Liberal PSI (group)","Core x liberal PSI (group)","Relative size","non-EPR group","Recent contestation","CPI","Polity (normalized)","GDP p.c. (logged)","Ethnic fractionalization","Recent contestation (country)","Recent contestation (TEK)","Only one non-core group","Age","Female","High education","Political interest"))


#######3.5 Independent variable alterations: power-sharing stocks instead of current levels (appendix 4.5, table A17)
#a) unsat_gov
m1_unsat_gov_psiage <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main1_stocks,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
m2_unsat_gov_psiage <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main2_stocks,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
#b) disc_grp
m1_disc_grp_psiage <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main1_stocks,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50))
m2_disc_grp_psiage <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main2_stocks,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50))
stargazer(m1_unsat_gov_psiage, m2_unsat_gov_psiage, m1_disc_grp_psiage, m2_disc_grp_psiage, type="text", style = "ajps", order=paste0("^", vars.order_stocks , "$"), omit = c("survey","region"), dep.var.labels = c("Government dissatisfaction", "Feeling discriminated"), notes = c("* p<0.1; ** p<0.05; *** p<0.01. Survey- and region-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Core","Included stock","Corporate PSI stock (group)","Core x corporate PSI stock (group)","Liberal PSI stock (group)","Core x liberal PSI stock (group)","Relative size","non-EPR group","Recent contestation","CPI","Polity (normalized)","GDP p.c. (logged)","Ethnic fractionalization","Age","Female","High education","Political interest"))


#######3.6 Individual moderation: political interest (appendix 4.6, table A18)#######
#a) unsat_gov
m3_unsat_gov_mod_interest <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main3_int,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
m4_unsat_gov_mod_interest <- glmer(as.formula(paste("unsat_gov", paste(c(re, re_i, re_i2,"survey",main4_int,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_unsat == 1 & group_s_number_unsat >= 50))
stargazer(m2_unsat_gov_mod_interest, m3_unsat_gov_mod_interest, m4_unsat_gov_mod_interest, type="text", style = "ajps")
#b) disc_grp
m3_disc_grp_mod_interest <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main3_int,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50))
m4_disc_grp_mod_interest <- glmer(as.formula(paste("disc_grp", paste(c(re, re_i, re_i2,"survey",main4_int,controls_gc, controls_i), collapse = " + "), sep = " ~ ")), family = binomial, nAGQ = 0, data = subset(grievances_ind, sample_disc == 1 & group_s_number_disc >= 50))
stargazer(m3_unsat_gov_mod_interest, m4_unsat_gov_mod_interest, m3_disc_grp_mod_interest, m4_disc_grp_mod_interest, type="text", style = "ajps", order=paste0("^", vars.order , "$"), omit = c("survey","region"), dep.var.labels = c("Government dissatisfaction", "Feeling discriminated"), notes = c("* p<0.1; ** p<0.05; *** p<0.01. Survey- and region-fixed effects included but not reported."), notes.append=FALSE, covariate.labels = c("Core","Corporate PSI (group)","Corporate PSI (difference, domestic)","Corporate PSI (difference, TEK)","Core x corporate PSI (group)","Core x corporate PSI (difference, domestic)","Core x corporate PSI (difference, TEK)","Liberal PSI (group)","Core x liberal PSI (group)","Relative size","non-EPR group","Recent contestation","CPI","Polity (normalized)","GDP p.c. (logged)","Ethnic fractionalization","Age","Female","High education","Political interest","Core x political interest","Corporate PSI (difference, domestic) x political interest","Core x corporate PSI (difference, domestic) x political interest","Corporate PSI (difference, TEK) x political interest","Core x corporate PSI (difference, TEK) x political interest"))



